//  Copyright (c) 2012 M.A. (Thijs) van den Berg, http://sitmo.com/
//
//  Use, modification and distribution are subject to the MIT Software License. 
// (See accompanying file LICENSE.txt or copy at http://www.stdfin.org/LICENSE.txt)

#ifndef STDFIN_DISTRIBUTIONS_GBM_FIRST_HITTING_TIME_HPP
#define STDFIN_DISTRIBUTIONS_GBM_FIRST_HITTING_TIME_HPP

#include <cstdlib> // abs
#include <cmath>
#include <limits>
#include <boost/math/distributions/complement.hpp>
#include <boost/math/distributions/normal.hpp>
#include <boost/math/constants/constants.hpp>

namespace stdfin {

template <class RealType>
inline RealType gbm_first_hitting_time_pdf(
        const RealType& S, 
        const RealType& mu, 
        const RealType& sigma, 
        const RealType& B, 
        const RealType& t) 
{
    RealType result;
    if (B==S) return std::numeric_limits<RealType>::infinity();
    if (t==0) return 0;
    
    RealType nu = mu / ( sigma * sigma ) - 0.5;
    result = std::abs( std::log(B/S) );
    result /= sigma * boost::math::constants::root_two_pi<RealType>() * std::pow(t, 1.5);
    result *= std::pow( B/S, nu);
    result *= std::exp( -0.5*nu*nu*sigma*sigma*t - std::pow(std::log(B/S),2) / (2.0*sigma*sigma*t) );
    
    return result;
}

template <class RealType>
inline RealType gbm_first_hitting_time_cdf(
        const RealType& S, 
        const RealType& mu, 
        const RealType& sigma, 
        const RealType& B, 
        const RealType& t) 

{
    RealType result;

    if (B==S) return 0; // direct hit
    if (t==0.) return 1; // because at this point B!=S)
    RealType eta = (B>S) ? -1 : 1;
    RealType lambda = std::abs(mu);
    RealType sigma_sqrt_t = sigma*std::sqrt(t);
    
    RealType z = std::log(B/S)/sigma_sqrt_t;
    RealType kappa = (mu - sigma*sigma/2.) / (sigma*sigma);
    
    // 1 - (C1*N1 + C2*N2) 
    // == 1 - C1 - C2 + C1*(1-N1) + C2*(1-N2) 
    RealType C1 = std::pow(B/S, kappa+lambda);
    RealType N1c = boost::math::cdf( boost::math::complement(boost::math::normal(), eta*(z+sigma_sqrt_t)) );
    RealType C2 = std::pow(B/S, kappa-lambda);
    RealType N2c = boost::math::cdf( boost::math::complement(boost::math::normal(), eta*(z-sigma_sqrt_t)) );

    result = 1 - C1 - C2  + C1*N1c + C2*N2c;
    return result;
}

} // namespaces
#endif
